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Abstract 


2 An ensemble Kalman filter (EnKF) is used in a suite of synthetic experiments to as- 

3 similate coarse-scale (25 km) snow water equivalent (SWE) observations (typical of satellite 

4 retrievals) into fine-scale (1 km) model simulations. Coarse-scale observations are assimi- 

5 lated directly using an observation operator for mapping between the coarse and fine scales 

6 or, alternatively, after disaggregation (re-gridding) to the fine-scale model resolution prior 

7 to data assimilation. In either case observations are assimilated either simultaneously or 
b independently for each location. Results indicate that assimilating disaggregated fine-scale 
s observations independently (method 1D-F1) is less efficient than assimilating a collection of 
io neighboring disaggregated observations (method 3D-Fm). Direct assimilation of coarse-scale 
a observations is superior to a priori disaggregation. Independent assimilation of individual 

12 coarse-scale observations (method 3D-C1) can bring the overall mean analyzed field close 

13 to the truth, but does not necessarily improve estimates of the fine-scale structure. There 
u is a clear benefit to simultaneously assimilating multiple coarse-scale observations (method 
is 3D-Cm) even as the entire domain is observed, indicating that underlying spatial error 
is correlations can be exploited to improve SWE estimates. Method 3D-Cm avoids artifical 
17 transitions at the coarse observation pixel boundaries and can reduce the RMSE by 60% 
is when compared to the open loop in this study. 

19 

20 Key words: assimilation, Kalman filter, downscaling, snow water equivalent, spatial cor- 

21 relation, multiscale 

22 
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1 1. Introduction 

2 Land surface data assimilation is mainly focused on surface temperature, soil moisture 

3 and snow. These variables interact with the atmosphere, which explains their direct impact 

4 on weather and climate predictions (Koster et al. 2004; Dirmeyer 2000). Land surface vari- 

5 ables have a large spatial variability that cannot be captured by the existing operational 

6 observing systems alone. Land surface models could help in either downscaling the available 

7 coarse satellite observations or interpolating the scattered point observations. Snow differs 
b from the other land surface states in its discontinuous local presence or absence, its cumu- 
s lative character and often long temporal autocorrelation length (Slater and Clark 2006). 

io For snow, a variety of observations can be assimilated. Firstly, there are numerous 
u snow water equivalent (SWE) or snowdepth point measurements, such as long term records 

12 from the U.S. Department of Agriculture’s (USDA) Natural Resources Conservation Ser- 

13 vice (NRCS) SNOwpack TELemetry (SNOTEL) network of snow pillow sites, the National 

1 4 Weather Service (NWS) Cooperative (COOP) weather stations network and short term, but 
is spatially dense, ground measurements from intensive field campaigns like the Cold Land 
is Processes Experiments (CLPX). Direct point snow data assimilation has been demonstrated 
17 by, e.g., Huang and Cressie (1996), Slater and Clark (2006) and Liston and Hiemstra (2008). 
is Point data have been widely used as validation for satellite data assimilation. 

19 Satellite observations of snow cover area or fraction and SWE-related quantities provide a 

20 second type of (partial) information about the snow state. As examples, the Advanced Very 

21 High Resolution Radiometer (AVHRR), Landsat Thematic Mapper (TM) and Moderate 

22 Resolution Imaging Spectroradiometer (MODIS) provide snow cover and albedo products 
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1 at fine resolutions (<1 km), but these visible or near infrared data cannot be collected 

2 at night or in cloudy conditions. Passive microwave sensor data do not suffer from these 

3 shortcomings, but have a much coarser resolution and come with large errors (Kelly et al. 

4 2003; Pulliainen et al. 1999). The Advanced Microwave Scanning Radiometer for the Earth 

5 Observing System (AMSR-E), the Scanning Multichannel Microwave Radiometer (SMMR) 

6 and the Spatial Sensor Microwace Imager (SSM/I) provide data at a 25 km resolution. Use 

7 of remote sensing data to update the SWE has focused on direct assimilation of MODIS snow 
b cover area data (Rodell and Houser (2004); Andreadis and Lettenmaier (2006); Zaitchik and 
s Rodell (2009)) or directly inverted (real or synthetic) SWE data from passive microwave 
io sensors (Sun et al. 2004; Andreadis and Lettenmaier 2006; Dong et al. 2007). 

u Thirdly, different approaches have been developed to merge complementary sources of 

12 remotely sensed snow data (Durand et al. 2008), to combine remote sensing products with 

13 ground measurements (Pulliainen 2006; Tart et al. 2000), and to assimilate a priori merged 

1 4 snow products into models (Liston et al. 1999). The merger of different observational sources 
is at multiple scales can also be achieved within the data assimilation framework, for exam- 

1 6 pie through the direct assimilation of multi-channel passive microwave brightness and near 

17 infrared data (Durand and Margulis 2006, 2007, 2008). 

is With few exceptions, notably Durand and Margulis (2007, 2008), snow data assimila- 

19 tion efforts have been limited to ID filtering, i.e. the observations and model units were at 

20 or brought to the scale of the model grid and assimilated independently for each grid cell. 

21 The integration of coarse-scale observations for land surface state estimation in a dynami- 

22 cal framework is often complex, because satellite data and computational model units (grid 

23 cells, catchments units, etc.) typically differ in scale of support. With more advanced as- 
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1 similation and scaling techniques, however, there is potential to extract information through 

2 downscaling satellite-scale observations into higher-resolution land model integrations. 

3 Coarse observations can be re-gridded to the finer modeling units prior to data assimi- 

4 lation. Examples include retrieving SMMR SWE over catchment units (Dong et al. 2007), 

5 and reconstructing MODIS snow cover and AMSR-E SWE to match the modeling scale 

6 (Andreadis and Lettenmaier 2006). Additionally, many synthetic studies take this approach, 

7 such as Sun et al. (2004) for SWE assimilation. More sophisticated a priori disaggregation 
b approaches based on statistical relationships have been used for soil moisture (Dubayah et al. 
s 1997; Kumar 1999; Crow and Wood 2002; Parada and Liang 2003b; Merlin et al. 2006). In 
io any case, the disaggregated products could then be useful for dynamical assimilation into a 
a finer-scale land model. 

12 Alternatively, the scale discrepancy can be dealt with more effectively and dynamically 

13 within the data assimilation framework by using a scaling observation operator. The dis- 

14 aggregation of coarse-scale information is then based on spatial error correlations that are 
is modeled within the assimilation system. For example in land surface modeling, Reichle 
is et al. (2001), Caparrini et al. (2004), Durand and Margulis (2007) and Zaitchik et al. (2008) 
17 used an averaging of a number of fine-scale forecasts to generate observation predictions at 
is the coarse observation scale, respectively for assimilation of soil moisture (synthetic L-band 

19 passive microwave), land surface temperature (AVHRR, SSM/I, Geostationary Operations 

20 Environmental Satellite GOES), SWE (synthetical AMSR-E) and water storage (Gravity 

21 Recovery and Climate Experiment, GRACE). 

22 In this paper, techniques for downscaling coarse-scale SWE observations to the under- 

23 lying fine-scale state variables within data assimilation is described and tested for the first 
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1 time within the National Aeronautics and Space Administration (NASA) Land Information 

2 System (Kumar et al. 2006; Peters-Lidard et al. 2007; Kumar et al. 2008a, b), version 5.0 

3 (LIS5.0). One of the challenges in assimilating coarse observations over a large domain of 

4 fine-scale model units is in the ensemble estimation of the forecast error covariance matrix, 

5 where there is often a problem with spurious long range correlations. Therefore, ensemble 

6 filter techniques mostly include some type of localization (Keppenne and Rienecker 2002; 

7 Ott et al. 2004; Hunt et al. 2006; Reichle and Koster 2003), which will be discussed in this 
b paper as well. As an alternative, multiscale Kalman filters could provide a way to replace 
s the forecast error covariance matrix with a multiscale tree (Zhou et al. 2008; Parada and 
io Liang 2004, 2008; Pan et al. 2009). 

a A number of ensemble Kalman filter (EnKF) approaches are explored that combine fine- 

12 scale land surface model simulations and synthetic satellite-scale observations. The focus 

13 will be on (i) the difference between ID (point) and 3D (spatial) filtering and (ii) the study 

14 of different 3D filter approaches in view of subpixel snow variability estimation. Reichle and 

15 Koster (2003) demonstrated that the 3D EnKF was superior to ID EnKF in a synthetical 
is soil moisture estimation study. In that study, observations and modeling units corresponded 
17 in scale of support (catchments), and the advantage in the 3D filtering was primarily in 
is data-sparse regions, because of the horizontal information propagation through the spatial 

19 forecast error structure (Hamill and Snyder 2000). This advantage of 3D filtering is also 

20 illustrated by Houser et al. (1998) and De Lannoy et al. (2009). In the synthetic experiment 

21 discussed here, however, observations are available for the entire domain and the focus is on 

22 the spatial disaggregation of coarse-scale observational information. 
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i 2. Experiment Setup 


2 Different filters are tested in an identical twin experiment. Snow pack evolution is sim- 

3 ulated at the fine-scale (1 km) as a reference “truth” for validation. A second integration 

4 with degraded forcing data mimics forecast modeling errors. Synthetic observations are 

5 generated from the truth integration at the coarse satellite-scale (25 km) by averaging the 

6 fine-scale snow pack and adding observation error. The coarse observation assimilation into 

7 the finer-scale degraded model simulations is then validated against the fine-scale truth. 

b The study domain includes the North Park area and part of the Rabbit Ears area in 
s Colorado, USA, which have been intensively monitored during CLPX-I. The rectangular 
io study area (left bottom corner: 40.255 N, -106.745 W; upper right corner: 40.995 N, -105.755 
a W) includes 75 xlOO fine-scale grid cells (0.01°~1 km) covering a central valley surrounded 

12 by mountains (Figure 1). 

13 SWE is simulated with the land surface model of the National Centers for Environmental 

1 4 Prediction (NCEP)/Oregon State University/ Air Force/Hydrologic Research Lab (Noah, 
is version 2.7.1, Ek et al. (2003)) at a 0.01° (~1 km) resolution for the winter of 2002-2003, 
is after a spin up starting in 2000. The NASA LIS5.0 (Kumar et al. 2008b) is used and 
17 adapted for this study. The forcings for the truth run are obtained from the North American 
is Land Data Assimilation System (NLDAS, original resolution 1/8°, Cosgrove et al. (1999)). 

19 For all simulations in this paper a fixed temperature lapse rate of -6.5 K.km -1 is used 

20 to disaggregate temperature data (Mitchell et al. 2004). The LSM parameters include a 

21 MODIS-based landmask and land cover parameters (1 km), the United Sates Geological 

22 Survey’s (USGS) GTOPO30 elevation map (~0.01°) to correct NLDAS forcings at the 1 km 
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scale, the Penn State University-USDA State Soil Geographic Database (STATSGO, 1/120°), 

2 a snow-free albedo product and a Noah-specific maximum snow albedo product from NCEP 

3 at 1 °, National Oceanic and Atmospheric Administration (NOAA)/AVHRR-based greeness 

4 fraction (0.144°), and an NLD AS-based climatological bottom temperature (1/8°). The 

5 default Noah soil, vegetation and general parameter tables are used. 

6 Synthetic SWE observations that mimic AMSR-E SWE retrievals are derived by aggre- 

7 gating the 1 km reference truth SWE simulation into twelve coarse-scale (25 km x 25 km) 
b grid cells (3 rows, 4 columns, see Figure 1). Random noise is added to mimic observation 
s errors and drawn from a normal distribution with a SWE-dependent standard deviation. 
io Specifically, we use cr 0 b s — 5 + 0.095 • SWE [mm]. The perturbation standard deviation is 
u minimally 5 mm, when no SWE is observed, and linearly increases to 100 mm for 1000 mm 

12 SWE. The average observation error magnitude is in the range reported for real passive 

13 microwave data, i.e. 5 to 45 mm in non-forested areas, increased by 5-10 mm over forests 

1 4 (Pulliainen et al. 1999; Derksen et al. 2003). Obviously, the assumed error function is only 
is a crude approximation of actual errors. In practice, the type of snow emission model (Kelly 
is et al. 2003; Pulliainen et al. 1999), local conditions and a number of physical details (includ- 
17 ing grain size, presence of liquid water, vegetation, Foster et al. (2005); Dong et al. (2005)) 
is will strongly influence the observation error in satellite-based passive microwave SWE prod- 

19 ucts. Furthermore, the AMSR-E products are known to have limited sensitivity to SWE 

20 for both very thin and deep snow packs (Dong et al. 2005). Nevertheless, we keep the 

21 whole range of synthetic observations, to illustrate the potential of the different assimilation 

22 techniques. 

23 A second set of Noah model forecasts (open loop) are generated using forcings from the 
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1 Global Data Assimilation System (GDAS), which is originally at a Gaussian grid (T170~0.7° 

2 until 29 October 2002, T254~0.5° until 31 May 2005). The forcings are of major impor- 

3 tance for the correct characterization of the SWE evolution (Slater et al. 2001; Mote et al. 

4 2003), and their errors will probably contribute most to SWE forecast errors in real experi- 

5 ments. Because the GDAS forcings result in heavily underestimated SWE simulations over 

6 the selected study domain, the precipitation is multiplied by a factor 3 to correct for the 

7 otherwise large forecast bias during the assimilation. This correction factor was determined 
b through trial and error. More sophisticated bias correction methods such as matching of the 
s cumulative distribution functions (Reichle and Koster 2004) or dynamic bias estimation (De 
io Lannoy et al. 2007) could be used, but are beyond the scope of this paper, which focuses on 
a the estimation of fine-scale spatial variability. The precipitation, air temperature, shortwave 

12 and longwave radiation, as well as the state variables are perturbed to generate ensembles 

13 for the assimilation, as will be discussed below. 

14 The satellite-scale (25 km) synthetic observations are assimilated into the fine-scale 
is (1 km) degraded model simulations from 30 September 2002 to 30 June 2003 (winter 2002- 

16 2003), either twice per month (every 15th and 30/28th of the month), or daily, to study the 

17 impact of the assimilation frequency on the analyses. To mimic the AMSR-E overpasses, 
is the assimilation time is set to 8:00 am UTC (1:00 am local). The night overpass is selected 

19 (Durand and Margulis 2007), because mid-day SWE retrievals are typically affected by liquid 

20 (melted) water in the snowpack. 
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3. Data Assimilation 


2 a. Ensemble Kalman Filter 

3 The EnKF is a Monte-Carlo variant of the Kalman filter (Evensen 2003). The idea behind 

4 the EnKF is that a small ensemble of model trajectories captures the relevant parts of the 

5 forecast error structure. Each member of the ensemble experiences perturbed instances of 

6 the observed forcing fields (representing errors in the forcing data) and is also subject to 

7 randomly generated noise that is added to the model parameters and prognostic variables 
b (representing errors in model physics and parameters). The error covariance matrices that 
s are required for the filter update can then be diagnosed from the spread of the ensemble 
io at the update time. The EnKF is flexible in its treatment of errors in model dynamics 
u and parameters. It is also very suitable for modestly nonlinear problems and has become 
12 a popular choice for land data assimilation (Andreadis and Lettenmaier 2006; Durand and 
u Margulis 2008; Kumar et al. 2008b; Pan and Wood 2006; Pauwels and De Lannoy 2006; 
i 4 Reichle et al. 2002a, b; Zhou et al. 2008). 

is The EnKF works sequentially by performing in turn a model forecast and a filter update, 
is To reflect the uncertainty in the state forecasts, perturbations of the model forcings and the 
17 initial state estimate are applied at each time step to generate N ensemble forecasts xj~ 
is (j — 1, • • • , N) with the model f: 

= f(x£ 1 ,u i ,w f_i), (1) 

19 where i denotes time, is the analysis from the previous time step (see below), u* repre- 

20 sents the forcings, and denotes the perturbations to the j-th ensemble member. The 
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1 model forecast x“ is given by the ensemble mean. As a convention, we use bold upper case 

2 symbols to refer to two-dimensional matrices, bold lower case symbols for one-dimensional 

3 vectors and non-bold italic symbols for scalars. 

4 When observations y * are available at time i, each ensemble member j is updated indi- 

5 vidually to obtain the a posteriori state estimate (or analysis): 

M + - M~ + K i[yi - h i(*ni- ( 2 ) 

6 where y{ is a suitably perturbed observation vector. The analysis x+ is again given by 

7 the ensemble mean. The function h* is the observation operator that maps the state to 

8 the observation space, and hj(x?“) denotes the model’s observation prediction for ensemble 

9 member j. The Kalman gain K* is identical for all ensemble members and determined 

10 by the (sample) error covariance Cov[x“, hj(x“)] between the forecast and the observation 
u predictions , the (sample) error covariance Cov[hj(x“), hj(x“)] of the observation predictions, 

12 and the observation error covariance Rji 

K = Cov[x“ , h i(k ~ )] [Cov[hj(x“), h^x" )] + R*] _1 . ( 3 ) 

13 For a linear (or linearized) observation operator Hj, i.e. hj(x“) ~ Hjic“, the sample error 

14 covariances can be written as: 

Cov[x“ , hi(x“)] = P, - Hf and Cov[h j (x-),h i (x")] = H ( P,-Hf, ( 4 ) 

15 where P“ is the forecast error covariance and superscript T denotes the matrix transpose, 
is In this paper we use the EnKF modules of the NASA Global Modeling and Assimilation 
17 Office (GMAO, Reichle et al. ( 2009 )) within LIS 5.0 (after modifying the LIS 5.0 implemen- 
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1 tation to enable the GMAO EnKF capability for horizontal error correlations, covariance 

2 localization, and distributed filtering). In each fine-scale (1 km) grid cell, the Noah model 

3 simulates a single snow layer with two prognostic variables, one for total column snow water 

4 equivalent (swe) and one for snow depth ( snd ). The EnKF is used to update these two 

5 model prognostic variables at the 1 km scale. Dropping the time index i, equation ( 2 ) and 

6 ( 3 ) can then be re-written for a given fine-scale grid cell k as 



K 


( swei 
\snd{ 




+ Kfy- 7 — y j ~] 


Cov ife) ’*"] [Cov[y-,y1 +R] _1 


( 5 ) 

( 6 ) 


7 with y- 7 perturbed SWE observations that are used to update swe J k ~ and snd J k ~, y- 7- = h (x- 7- ) 

8 is short-hand for the corresponding model predictions of these observations, and R denotes 

9 the corresponding observation error covariance. 

10 The first term in the expression for the Kalman gain includes the error correlation between 
a forecasted snow depth ( snd ~ ) and SWE (y“). This correlation information is the basis for 

12 updating snow depth ( snd ~ ) in response to SWE observations (y). If to observations are 

13 used to update SWE and snow depth for a given fine-scale grid cell, the Kalman gain K in 

14 equation ( 5 ) is a 2 x to matrix, and the observations y and observation predictions y - are 
is column vectors with to elements. In sections c and d, we will specify how exactly the coarse- 

16 scale SWE observations are used (to update the fine-scale model prognostic variables swe k 

17 and sndjT) by relating y- 7 and y- 7- to observation and model variables. In any case, successful 
is downscaling of coarse-scale observational information relies on ancillary information at the 
i9 smaller scale, for example from micrometeorological or terrain information, and on horizontal 
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1 error correlations that enable horizontally distributed updates. 

2 b. Ensemble Perturbation 

3 Random perturbation fields (normally distributed, zero mean, spatially correlated, tem- 

4 porally uncorrelated) are added hourly to the longwave radiation (LW, standard deviation, 

5 stdv — 10 W.m -2 ) and near surface air temperature (AT, stdv — 1 K) forcings, as well 

6 as to the forecasted SWE ( stdv — 0.0025 m) and snowdepth ( stdv — 0.01 m). The pre- 

7 cipitation (P) and shortwave radiation (SW) are perturbed through multiplication with a 
b random lognormally distributed field with stdv — 0.2 and stdv — 0.5, respectively, for the 
s standard deviation of the multiplication factor. Cross-correlation between the forcing per- 
io turbations is included ( LW-SW : -0.3, SW-P: -0.1, SW-AT : 0.3, LW-P: 0.5, LW-AT : 
u 0.6, P-AT: -0.1), as well as cross-correlation between swe~ and snd~ forecast state (0.9) 

12 perturbations. The magnitude of the perturbations is chosen in a realistic range and checked 

13 to allow near-optimal filter performance (see later). 

14 Experiments are conducted with several different spatial correlation lengths (L corr) of 
is the perturbation fields. Values of L corr range between approximately 5 km and 100 km (i.e., 
is 0.05°, 0.1°, 0.2°, 0.3°, 0.5°, 0.7°, 1°) and are always chosen to be identical for all different 
17 perturbation fields in a given experiment. When the spatial correlations are not included 
is in the Kalman gain calculation, then the choice of L corr does not matter and the spatial 

19 correlations are effectively forced to zero in the Kalman gain calculation. The number of 

20 ensembles is set to N — 12. 
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1 c. Observation Pre-processing: A Priori Observation Disaggregation (1D-F1, 3D-Fm) 

2 Perhaps the simplest approach for the assimilation of coarse-scale observations into finer- 

3 scale model simulations is to create a surrogate observational grid with the same resolution 

4 as the fine-scale model grid. Here, the simplest possible disaggregation operation is used by 

5 assigning to each fine-scale (1 km) grid cell the value of the observation in the coarse (25 km) 

6 grid cell that contains the fine-scale grid cell. 

7 An easy assimilation approach is then to apply the EnKF independently at each fine- 
b scale grid cell k, using a single disaggregated and collocated observation swe°^ s as illustrated 
s in Figure 2 . In this approach, spatial forecast error correlations are entirely disregarded. 
io This filter is referred to as 1 D-F 1 , a one-dimensional (ID) filter that uses exactly one fine- 
u scale observation (FI) for the analysis increment computation at a given fine-scale grid cell. 

12 Formally, we compute the analysis increments K[y J — y- 7- ] for the entire domain by looping 

13 through all fine-scale grid cells k and using (for each k ) 

y j = swef 8 ' 3 ( 7 ) 

y j ~ = swe? k ~ (8) 

14 as inputs to equation ( 5 ). 

15 A straightforward extension of the 1 D-F 1 approach is to use several (to) neighboring 
is observations (after disaggregation to the fine scale) in the analysis increment computation 
17 for a given fine-scale grid cell k (Figure 2 ). This approach limits the edge-effect at the 
is transition line between neighboring coarse observation pixels through the consideration of 
i9 horizontal error correlations in the Kalman gain computation. We refer to this approach as 
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1 3D-Fm, because horizontal error correlations are taken into account in a “three-dimensional” 

2 (3D) filtering approach (following the convention of Reichle and Koster (2003)) and because 

3 m (disaggregated) fine-scale (Fm) observations are used to update a given fine-scale grid 

4 cell. Formally, we compute the analysis increments for the entire domain by looping through 

5 all fine-scale grid cells k and using (for each k ) 


y j 



SWe °k2 ,j 

\swefy) 

/Wr\ 

swe U 

\ swe ij 


( 9 ) 


( 10 ) 


e as inputs to equation (5), where k 1: k-z, ..., k m indicate the m (fine-scale) observations and 

7 predictions that are included in the update of fine-scale grid cell k. The choice of the 

8 observation selection area (influence radius around the analysis grid cell k ) determines the 
s observation dimension m. Because a larger influence radius also means a more complicated 
io inversion for the Kalman gain calculation, we limit the influence radius to 10 km when 
u 1 km resolution observations are used (maximum m — 317, m is smaller near the domain 

12 boundaries). The observation error covariance matrix R* is assumed diagonal. Obviously, 

13 as the influence radius approaches zero, the 3D-Fm approach reduces to the 1D-F1 method. 
M In our synthetic experiment, the dynamic observation error variance is known for the 
is coarse-scale observations, but not necessarily for the disaggregated observations. The ob- 
is servation error includes both the instrument error and the representativeness error of an 
17 observation operator. While the first error term is typically assumed to be white Gaussian, 
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1 the second one can affect the actual observation error magnitude and correlations for the 

2 disaggregated observations. For simplicity, we do not consider spatially correlated observa- 

3 tion errors and the same observation error magnitude is used for both fine- and coarse-scale 

4 observation assimilation. 

5 d. Observation Operator: Forecast Upscaling (3D-C1, 3D-Cm) 

6 It is more appealing to avoid the observation disaggregation prior to assimilation and 

7 to perform the downscaling within the filtering algorithm instead (Reichle et al. 2001; Du- 

8 rand and Margulis 2007; Zaitchik et al. 2008). Two disparate grids are then present in the 
s multiscale filtering algorithm: the coarse observation grid and the fine-scale simulation grid, 
io indicated by Greek and Roman letters, respectively. For example, a coarse-scale observation 
a is denoted with swe^ 8 . 

12 A first option is to use only the single coarse-scale observation swe°£ 3 that covers the fine- 

13 scale grid cell k to update swejT and snd J k ~ with equation (5) as illustrated in Figure 2. In 

14 this case, the coarse-scale observation predictions are computed as the appropriate average 

15 over the corresponding fine-scale model predictions. Formally, we have 

y- 7 = swe (11) 
y>- = 1/625 Egwoef ( 12 ) 

1 6 where ki , Z = 1,2, ...,625 indexes the 25 x 25 fine-scale (1 km) grid cells contained within 

17 the coarse-scale (25 km) grid cell k. This approach is referred to as the 3D-C1 EnKF 
is here, because it uses one coarse-scale observation pixel to update each fine-scale state. The 
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1 error correlation between the fine-scale model state (to be updated) and the coarse-scale 

2 observation prediction (first term of the Kalman gain in equation (5)) then downscales the 

3 observational information into fine-scale updates. 

4 The 3D-C1 technique can be expanded to include surrounding coarse observations in 

5 the analysis, thereby using spatial forecast error correlations that overarch the boundaries 

6 between the coarse observation pixel areas. Now, multiple (to) coarse-scale observations are 

7 used. Formally, equation (5) will be used with 


y j 



/swe°£ s ’ j \ 

swe°£ s ’ j 

WzS'V 

1/625 Egotoe^X 
1/625 Egswei- 

/ 1/625 EJ l\sw4jj 


(13) 


(14) 


8 where K m , m — 1,2, ...,12 denotes a coarse-scale grid cell in our study domain and ki itn , 

9 l — 1,2, ...,625 indexes the fine-scale grid cells contained within K m . This technique is 

10 referred to as 3D-Cm, involving multiple coarse-scale observations and in a 3D filtering 
u approach. As in Reichle and Koster (2003) we suppress spurious long range correlations 

12 in the forecast ensemble through (element-wise) Hadamard multiplication of the sample 

13 error covariance terms in equation (5) with a distance-dependent and compactly-supported 
M function. Specifically, we use the 5th-order polynomial equation (4.10) of Gaspari and Cohn 
is (1999) and set the compact support scale equal to 2.5 times the spatial correlation length 
is l-corr that is used for the ensemble spatial perturbation fields (i.e. 2.5°, 1.75°, 1.25°, 0.75°, 
17 0.5° and 0.25°). In practice, the localization means that observations beyond a given distance 
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1 from a specific fine-scale grid cell cannot affect the analysis at that grid cell. 

2 Equation (5) is written for a single fine-scale grid cell k but involves knowing a potentially 

3 large number of fine-scale model states for the observation prediction ( y J “) computation. 

4 For implementation on a parallel architecture, parallelizing the loop through k is thus not 

5 necessarily the best strategy. Here, we implemented the 3D-C1 and 3D-Cm update methods 

6 by looping through coarse-scale grid cells and simultaneously updating all 625 fine-scale grid 

7 cells within a given coarse grid cell. 

b 4. Results 

s The assimilation results are validated against the reference truth at the fine scale. Unless 

10 otherwise specified, the results are analyzed for assimilation twice a month (19 events be- 

11 tween 30 September 2002 and 30 june 2003) over the entire domain of 7500 fine-scale 1 km 2 

12 simulation grid cells (covered by 12 coarse-scale observation pixels). 

13 a. SWE Open Loop Integration and Coarse Observations 

14 We first examine the synthetic coarse-scale observations and the (ensemble mean) open 
is loop (no assimilation) integration. The top three rows of Figure 3 show five snapshots of 
is the SWE field for the truth integration, the observations, and the open loop integration. 
17 Figure 4 shows the root-mean-square error (R.MSE) of SWE (in time and in space) along 
is with a time series of the domain average SWE. By construction the domain average SWE 
19 for the observations very closely approximates the true domain average SWE. Note also that 
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1 the open loop fields for the various choices of perturbation correlation lengths (L corr) are 

2 nearly identical, thus we do not distinguish between them for plotting. The observations 

3 show a smaller RMSE than the open loop integration. Because the observations are based 

4 on the truth simulation, there is no bias in the observations. For the open loop integrations, 

5 by contrast, some local or temporal bias (in addition to random errors) cannot be avoided. 

6 Table 1 summarizes the space-time RMSE over the 19 analysis time steps. The RMSE 

7 for the observations is about 44 mm. As mentioned above, the ensemble mean open loop 
b simulations for different L corr perform similarly, with an RMSE of about 77 mm. 

s The difference between the truth and the open loop is caused by the difference in NLDAS 
io and (scaled) GDAS forcings. Even though the GDAS forcings are at a very coarse resolu- 
u tion, the modeled SWE distribution is quite reasonable. This is because the temperature 

12 correction (lapse rate of - 6.5 K.km -1 ) downscales the temperature allows a fine-scale dis- 

13 crimination between snowfall and rainfall. (But recall that a multiplication factor of 3 was 

14 used for the precipitation to bring the GDAS snowfall input to a more realistic level.) The 
is peak in the spatial RMSE for both the open loop estimates and observations occurs during 
is the melt period. This is caused by an increasing patchiness and variability in the actual 
17 snow amounts (as evidenced in the last column of Figure 3 ), which is impossible to reflect in 
is the coarse-scale observation and also not fully captured by the open loop integration. The 

19 open loop shows almost a month delay in melting off all snow during early spring (Figure 4 ), 

20 because of an accumulated snowpack error. 

21 Figure 5 shows a detail of the temporal SWE evolution for two individual fine-scale 

22 pixels. At the fine-scale, the coarse observations can differ significantly from the truth (e.g. 

23 during the snow accumulation in pixel ( 49 , 49 )), because they do not represent the local 
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1 fine-scale SWE evolution. This fine-scale observation bias is important to keep in mind for 

2 ID filtering with disaggregated observations. Note also that the SWE forecast uncertainty 

3 (P“ as measured by the ensemble spread; not shown) increases from zero when no snow is 

4 present to a maximum value during deep snow packs, and then drops back to zero again 

5 after snowmelt. 

6 b. SWE Analyses 

7 1) Space-Time Averaged Statistics 

b For most filter scenarios the space-time average RMSE (Table 1) is lower than for either 
s the open loop (77 mm) or the observations (44 mm). This indicates that both the spatial 

10 mean field and the fine-scale variability are improved and supports the premise that assim- 

11 ilation products are better than either forecasts or observations alone (Reichle and Koster 

12 2005). 

13 The 1D-F1 and 3D-Fm results show a similar RMSE (40-44 mm). Larger reductions 

14 in RMSE can be observed for the 3D-C1 and 3D-Cm with Lcorr > 20 km. For just 19 
is assimilation events, the 3D-Cm approach achieves an RMSE decrease up to 60% when com- 

1 6 pared to the open loop and up to 25% when compared to the observations. For 3D-C1, 

17 the RMSE values are slightly higher (around 38 mm) than for 3D-Cm (around 34 mm) for 
is Lcorr > 20 km, but for shorter Lcorr the 3D-Cm shows less assimilation impact. With daily 

19 assimilation, the RMSE over the same 19 time steps could not be significantly reduced. In 

20 the following subsections we analyze these findings in more detail. 
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1 2) Temporal and Spatial SWE Patterns 

2 Figure 3 shows the spatial analysis SWE pattern for the different filter approaches. The 

3 time series of spatial R.MSE for the assimilation analyses show a relatively smooth evolution 

4 (Figure 4). This means that the updates are accepted without significant drift in the spatial 

5 mean and that the snowpack has a considerable memory. During early snow accumulation, 

6 none of the filters, except the 3D-Fm (see dip in the spatial R.MSE evolution in Figure 4 and 

7 spatial field in Figure 3 on 15 October 2003), is able to remove the excessive snow to match 
b the true absence of snow. A very thin layer is kept, which is generated through the ensemble 
s perturbation. During this initial period with extremely low snow amounts the update is 
io negligible, because the model spread is still negligible, while the observation uncertainty has 
u a set minimum value of 5 mm. During the ablation period, the open loop snow melt delay 

12 can be largely reduced (Figure 4) as a result of earlier snow pack corrections and the few 

13 instantaneous updates during the melt. The melt delay cannot be reduced for the filter 

14 scenarios relying excessively on the model predictions and hence resulting in analyses close 
is to the open loop, like 3D-C1 with Lcorr — 0 km and 3 D-Cm with very low Lcorr and a small 
is localization scale. With a limited spatial error correlation, the error covariance between a 
17 single fine-scale grid state and the observation prediction (based on a number of fine-scale 
is grid elements) is low and the Kalman gain is limited. 

19 In the spatial dimension, a clear coarse block structure in the analyses is evident for 

20 the 1 D-F 1 , 3D-Fm and 3D-C1 filters (Figure 3). The 1D-F1 strongly forces all fine-scale 

21 state values to the observed coarse-scale SWE. Towards the melt season, the forecasts have a 

22 reduced spread and the spatial structure of the model simulations becomes apparent. For the 
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1 3D-Fm with an influence radius of 10 km, the coarse pixel boundaries are slightly smoothed, 

2 but they are still very obvious. By design, the 3D-C1 will always show the coarse pixel 

3 boundaries, because the state variables within a coarse area are all updated using only the 

4 single overlying coarse observation. However, when Lcorr— 0 km, the boundaries are less 

5 noticeable, because of a negligible observation impact (see above). By construction, the 

6 3D-Cm does not show coarse observation transitions in the analysis result, since it uses a 

7 collection of coarse observations centered on each individual analysis point. 

b The spatial pattern of the R.MSE (Figure 4) again shows the block patterns for some 
s filters. Here, it is also clear that the spatial variability in the R.MSE reduces with a longer 
io Lcorr for 3D-C1 and 3D-Cm, and fine-scale locations with particular high R.MSE values in 
u the open loop integration are generally much improved through assimilation. 

12 3 ) Spatial Forecast Error Correlation and Localization 

13 The results for the 3D-Cm and the 3D-C1 perform poorly when Lcorr < 20 km (Table 1). 

14 No drastic change in performance is found with Lcorr > 20—30 km. This threshold distance 
is is roughly equal to the dimension of the coarse-scale observation, and also approximates 
is the mean of the observed precipitation error correlation length (as measured by analyzing 
17 the spatial correlation length in the NLDAS minus GDAS precipitation difference fields), 
is although the latter was found to vary significantly in time (not shown). 

19 For 3D-C1, updates near the borders in the coarse observation pixel may be limited, 

20 when Lcorr is smaller than the pixel dimension. In the border areas, the average distance to 

21 all other fine-scale pixels within the coarse observation area is longer than for center pixels, 
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1 and the cross-correlation with fine-scale variables at longer distances is smaller, especially 

2 when the distance reaches beyond L corr. This causes a limited covariance Cov[x“, hj(x“)] 

3 between the forecast error in a single fine-grid location and an observation prediction that 

4 is based on a number of fine-scale grid forecasts. The sample observation prediction error 

5 covariance Cov[hj(x“), hj(x“)] (which is independent of the fine-scale update location) is 

6 also smaller for shorter L corr. Therefore, a smaller L corr results in a smaller gain factor, 

7 mainly for the fine-scale analysis closer to the borders. For L corr > 20 — 30 km, the forecast 
b error correlations are near-uniformly high over the each individual coarse satellite pixel area, 
s Increasing the L corr beyond the coarse pixel dimension adds little to the covariances in 
io the Kalman gain for each fine-scale update within the coarse pixel. This may explain why 
u Durand and Margulis (2008) also found a maximum assimilation efficiency for an exponential 

12 variogram correlation length of 25 km when assimilating 25 km resolution SWE-related 

13 observations, without significant degradation for longer L corr. 

u For the 3D-Cm approach, multiple coarse observations are included in each fine-scale 
is update and a covariance localization scale equal to 2.5 • L corr is applied to limit the impact 
is of more distant observations. For small L corr (relative to the domain size), this avoids 
17 adverse effects caused by spurious long range correlations. However, the smaller L corr plus 
is the localization limit the magnitude of the covariances in the calculation of the Kalman gain. 

19 This results in analyses close to the open loop simulation (compare the spatial and temporal 

20 RMSE in Figure 4 and see the analysis R.MSE in Table 1 for L corr — 5 and 10 km). For 

21 longer Lcorr , the corresponding longer covariance localization scale only marginally reduces 

22 the sampled spatial correlations in this bounded study area. With longer Lcorr, the update 

23 is increased locally and spread to other locations as well (smoothing), because of larger 
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1 forecast ensemble-based error covariances in the calculation of the Kalman gain. 

2 In our study, a spatially isotropic error correlation and localization is imposed on the 

3 random fields, but the error field may not exactly represent the actual forecast error structure. 

4 Given the large elevation gradients and forecast errors introduced by different forcings, a 

5 complex error field can be expected. A good characterization of the error field and a region- 

6 dependent localization may further improve the results. 

7 4) Fine-Scale Subpixel SWE Variability 

b Figure 6 shows scatter plots of the fine-scale assimilation estimates (forecasts and anal- 
s yses) within a single coarse grid cell (2,2) versus the true values just before and just after 
io one particular assimilation update. The assimilated coarse-scale observation is also shown, 
a These plots explain how each assimilation algorithm has a different effect on the spatial 

12 mean value and spread in the fine-scale SWE analyses. With the 1D-F1 filter, all fine-scale 

13 SWE values are drawn towards the same observed coarse-scale value. The spatial variability 

14 in the analysis result is much smaller than in the forecast. The same holds for the 3D-Fm 
is (not shown), but the analysis has slightly more spatial variability, because of the interaction 
is between neighboring innovations at each point. With the 3D-C1, the cloud of fine-scale pixel 
17 values is moved towards the observation, based on the coarse-scale difference between the 
is observation and the forecast, but the forecasted spatial variability is generally maintained. 

19 The perhaps most interesting finding is that for the 3D-Cm, the spatial variability of the 

20 analysis is strongly improved and better mimics the truth, indicated by the analysis scatter 

21 points lined up around the 1-1 fine. Apparently, including neighboring observations and a 
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1 reasonable forecast error structure can help to improve the fine-scale SWE estimation over 

2 using only the observation at the analysis point of interest. 

3 We also examined the spatial correlation coefficient of SWE fields from assimilation ver- 

4 sus the reference truth. The correlation coefficient measures the skill of the assimilation 

5 estimates in terms of fine-scale structure. Figure 7 shows the time series mean of the cor- 

6 relation coefficient for each coarse-scale pixel area and for a variety of filter scenarios. By 

7 construction, the spatial correlation coefficient between the synthetic observations and the 
b truth is negligible, because the observations are purely randomly perturbed. The filters with 
s disaggregated observations (1D-F1, 3D-Fm) degrade the spatial patterns with respect to the 
io open loop. Overall, the 3D-C1 has only a limited positive effect and the spatial structure is 
a only better than the open loop when Lcorr > 20 — 30 km. However, the 3D-Cm scenarios 

12 with larger Lcorr show a great correlation improvement, mainly in the center coarse areas 

13 that benefit most from surrounding information and a relatively gentle topography (and 

14 hence relatively isotropic precipitation error correlation field). For coarse grid cell (2,2), for 
is example, the correlation coefficient increases from a mean of 0.15 for the open loop to 0.64 
is for the 3D-Cm with Lcorr — 20 km. Only in two grid cells ((3,2) and (3,3)) the 3D-Cm does 
17 not improve the spatial variability, which is most likely because of the simple forecast error 
is correlation structure in combination with a highly dynamic and spatially variable snow pack 

19 distribution in that area. 

20 Figure 5 shows the temporal SWE evolutions at two individual fine-scale points within 

21 the coarse region (2,2). Most of the analyses better approach the truth than either the 

22 open loop or the observations, but none of the ID and 3D filters are necessarily optimal at 

23 all locations within the coarse-scale pixel. At fine-scale grid cells where the disaggregated 


24 



1 coarse observations are consistently biased when compared to the fine-scale truth, 1D-F1 

2 assimilation will degrade the analysis results compared to the open loop (e.g. during the 

3 snow accumulation period in pixel (49,49)). The ID approach with disaggregated observation 

4 values has an observation bias problem at some individual locations. This problem can 

5 be overcome by assimilating the observations at the coarse-scale with a well defined error 

6 correlation structure and the proper scaling observation operator (as in 3D-C1 or 3D-Cm). 

7 For each update using multiple observations (3D-Fm, 3D-Cm), the SWE analysis for a given 
b fine-scale pixel can exceed (or be less than) the forecast and the local observation. This 
s is mainly obvious during the melt phase, when some patches still have a high SWE value 
io with a high forecast uncertainty, while other patches have a very shallow snowpack with a 
a very limited uncertainty. An innovation at the deeper snowpack in the influence area around 

12 an analysis point strongly impacts a location with a limited SWE (limited uncertainty and 

13 limited update by the collocated observation). In general, the 3D-Cm analyses are closest 

14 to the true reference simulation. For 3D-C1, most analyses are moved towards the same 
is direction with respect to the forecast, but that is not necessarily the best solution at each 
is fine-scale location. The 3D-C1 filter with small l-corr results in an analysis close to the open 
17 loop, because of the small Kalman gain, as discussed earlier. 

is 5) Filter Diagnostics 

19 Figure 8 shows for different assimilation cases the frequency distribution of the normalized 

20 ensemble mean innovations, i.e. each innovation value [y* — HjX^Jfc (at the fine scale) or 

21 [y* — HjX~] K (at the coarse scale) is normalized by the square root of its filter-estimated 


25 



1 


standard deviation |H,P-Hf + RilJ ' 2 or [H,P,-Hf + RJi ' 2 (fine or coarse scale, resp.). 

2 The normalized ensemble mean innovations should approximately obey a standard-normal 

3 distribution (Gaussian with mean zero and variance one), if the model is linear and the filter 

4 operates in accordance with its underlying assumptions. The innovation histograms provide 

5 a rough indication of whether the model and observation error parameters are appropriately 

6 chosen or not. The histograms show all innovations over space and time (for the 3D-C1 and 

7 3D-Cm, there are only 12-19 innovations, while for the 1D-F1 and 3D-Fm there are 7500-19). 
b The figure shows a small bias in the innovation distributions. This is caused by local biases 
s in the forcings. When a daily assimilation is performed during the period from 30 September 
io 2002 through 30 June 2003 (274 assimilation time steps instead of 19), the histograms (not 
a shown) exhibit similar features, but they show no bias. 

12 For most filter scenarios, the standard deviation of the histograms is close to 1, indicating 

13 that the filter parameters are near optimal. During early snow accumulation and late abla- 
M tion, the forecast uncertainty is bounded and it increases during deep snow packs. Similarly, 
is the observation error was increased with deeper snowpacks. This similar trend in observar 
is tion and forecast errors keeps the filter close to its optimal operation during the entire snow 
17 season. However, with real observations, it remains to be seen how the observation error 
is evolves in time and how it compares to the forecast uncertainty in each time period. Of all 

19 filter configurations, 3D-C1 with Lcorr — 0 km suggests the least optimal filter operation, 

20 with an excessive innovation spread. This is because the update is too limited to reduce the 

21 forecast uncertainty during filtering. Observation error variance sensitivity tests (not shown) 

22 indicate that scaling of the chosen dynamic observation error variance for the disaggregated 

23 observations could slightly reduce the R.MSE for the 3D-Fm filter (for the specified forecast 
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1 error). However, the corresponding innovation statistics show a suboptimal filter operation. 

2 This suggests that the better R.MSE estimates may have resulted from compensation for 

3 imperfectly defined model errors. 

4 Note that the innovation diagnostics are always available, including when satellite ob- 

5 servations are assimilated, and can be used to identify obviously inadequate assimilation 

6 parameters (such as Lcorr — 0 km for 3 D-C 1 ) in the absence of reliable validation data. 

7 5. Conclusions 

8 The assimilation of coarse-scale (25 km) water equivalent (SWE) observations into fine- 

9 scale (1 km) model simulations is studied through a number of synthetic experiments with 

10 GMAO EnKF variants implemented in LIS 5 . 0 . Both a 1 -dimensional (ID) or point filter 
u and a variety of 3 -dimensional ( 3 D) or spatial filters are explored in a synthetical study. The 

12 truth is generated at the fine scale by forcing the Noah model with NLDAS data. Synthetic 

13 observations are generated through aggregation of these fine-scale simulations and addition of 

14 observation error. To assure a reasonable coefficient of variation in the observation error and 
is in favor of the filter performance, the observation error variance is a function of the spatial 
is mean SWE amounts. A degraded open loop model integration is simulated by forcing the 
17 Noah model with coarse GDAS data and imposing different types of spatially correlated 
is random perturbations. 

19 It is shown that coarse satellite products can improve both the SWE spatial mean and 

20 variability estimation over the open loop simulations, or in other words, fine-scale model 

21 simulations can downscale coarse satellite products to extract useful information for the 
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1 fine-scale SWE estimation. Because of the temporal evolution of the SWE amounts over the 

2 winter, the forecast uncertainty is small during early accumulation and later during the melt- 

3 ing season. Therefore, during the first days of simulated snow accumulation, the very small 

4 simulation spread and the fixed minimum observation error limits the assimilation impact. 

5 The predicted month delay in snow melt can be mostly removed through assimilation. 

6 When the coarse observations are disaggregated prior to assimilation, a number of prob- 

7 lems arise. The ID filter pushes the analyses to the coarse observation value, which improves 
b the spatial mean SWE estimation but removes most of the subpixel spatial variability. Fur- 
s thermore, significant fine-scale observation bias can be observed. If a number of fine-scale 
io observations are used to update each analysis point in the 3 D-Fm, the spatial variability is 
a slightly better, but the computational cost becomes increasingly expensive in the inversion 

12 part of the Kalman gain calculation. For the actual assimilation of passive microwave prod- 

13 ucts (e.g. AMSR-E), disaggregating the coarse observations to finer scale observations may 

14 render the observations nearly useless, because the coarse products already have an unfa- 

15 vorable signal-to-noise ratio. The correct scaling of the observation error for disaggregated 
is observations needs further research to optimize the filter at the fine scale. 

17 The best approach is to assimilate the coarse-scale observations directly with a 3 D filter 
is and a properly defined forecast error correlation structure. For all 3 D filters using the coarse 

19 observations, it is found that the analyses are best when the spatial forecast error correlation 

20 length is equal to or larger than 20-30 km, which corresponds to the dimension of the coarse 

21 observation pixels, and also to the approximate correlation length of the precipitation error 

22 field. The results degrade significantly for shorter correlation lengths. When each fine-scale 

23 grid point is updated using the overlying coarse observation only ( 3 D-C 1 ), the spatial mean 
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1 SWE field can be improved over that of the observations or open loop simulations alone, 

2 and the spatial subpixel variability can be enhanced slightly. With additional inclusion of 

3 surrounding coarse observations (3D-Cm), there is a significant improvement in estimating 

4 the subpixel variability. Furthermore, artificial boundaries in the analysis field, caused by 

5 the boundaries in the coarse observations, are completely removed. 

6 Even though the entire study domain is observed, there is a substantial value in including 

7 coarse observations from neighboring areas to enhance the fine-scale SWE structure estimar 
b tion within the 3D-Cm filter. The SWE estimation can probably be further improved after 
s optimization of the spatial error structure and covariance localization for the forecast error 
io covariances in spatially complex terrain areas. In summary, the 3D-Cm filter is the method 
u of choice to assimilate coarse observations, because it (i) improves the spatial mean SWE 

12 analysis, (ii) substantially enhances the subpixel SWE variability estimation, (iii) avoids ar- 

13 tificial transitions at the coarse observation boundaries and (iv) is computationally no more 

1 4 expensive than any of the other filter approaches. 
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Fig. 1. Digital elevation model of the studied domain. The two plus-symbols mark locations 
for which individual time series are plotted in figure 5, with the upper one at fine-scale row 
and column (49,49) and the lower one at (27,46). 
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3D-Fm: multiple fine-scale obs 
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Fig. 2. Schematic of EnKF approaches illustrated for four coarse-scale pixels (gray shad- 
ing), each containing 7x7 fine-scale pixels. For the state update of a given fine-scale pixel 
(black filled square), we use observations within an influence area (white circle) and a corre- 
sponding set of fine-scale model forecasts (thick black fine). For 1D-F1 and 3D-Fm fine-scale 
(disaggregated) observations are used (small crosses). For 3D-C1 and 3D-Cm, coarse-scale 
observations are used (center location is indicated by a larger crosses). In the 3D- Cm exam- 
ple, the observation at the white cross location is outside the influence area and not used to 
update the marked fine-scale state variable. 
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Fig. 3. Snapshots of SWE fields for (columns) 15 Oct 02, 30 Nov 02, 15 Jan 03, 28 Feb 03 
and 15 Apr 03. (Top row) truth, (second row) synthetic observations, (third row) ensemble 
mean open loop forecasts and (rows 4-9) analyses obtained with several filter approaches. 
The spatial correlation length Lcorr in the forecast perturbations is indicated in parentheses. 
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Fig. 4. Temporal (field) and spatial R.MSE (bold black line) of the fine-scale analysis SWE 
for the same assimilation scenarios as in figure 3. All analysis and forecast time steps at 
8:00 UTC in the period between 30 September 2002 through 30 June 2003 are included. 
The dashed black line represents the total domain averaged SWE. The solid gray line in the 
upper plot for observations shows the true domain-average SWE. 
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Fig. 5. Observed SWE (Syn Obs), ensemble mean open loop (Ens OL) and analyses of 
different assimilation algorithms at two fine-scale locations (a, b) in the same coarse ob- 
servation area (see figure 1). The observations are identical for both fine-scale locations. 
Assimilation is performed twice each month. 
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Fig. 6. Comparison of the fine-scale (left) forecasts and (right) analyses with the true SWE 
distributions within a single 25 x 25 km 2 coarse pixel (second row, second column of coarse 
grid) at 1 particular time step (30 Apr 2003; 8:00 am UTC). The dashed fines indicate the 
spatial means and the box shows the assimilated coarse-scale observation value. 
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Fig. 7. Temporal mean of the spatial correlation coefficient between analyzed and true SWE 
as function of perturbations correlation length for (dashed) open loop, (black dotted) 3D- 
Cl, and (black solid) 3D-Cm. Methods 1D-F1 (Lcorr — 0 km) and 3D-Fm (Lcorr — 5 and 
10 km) are shown with + symbols. Each subplot corresponds to a coarse grid cell (Figure 1). 
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Table 1. R.MSE of SWE [mm] over 19 update times and 7500 fine scale locations. Lcorr 
is the spatial correlation length in the forecast perturbation fields. 


l-corr 
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An ensemble Kalman filter (EnKF) is used in a suite of synthetic experiments to 
assimilate coarse-scale (25 km) snow water equivalent (SWE) observations (typical of 
satellite retrievals) into fine-scale (1 km) model simulations. Coarse-scale observations 
are assimilated directly using an observation operator for mapping between the coarse 
and fine scales or, alternatively, after disaggregation (re-gridding) to the fine-scale model 
resolution prior to data assimilation. In either case observations are assimilated either 
simultaneously or independently for each location. Results indicate that assimilating 
disaggregated fine-scale observations independently (method 1D-F1) is less efficient than 
assimilating a collection of neighboring disaggregated observations (method 3D-Fm). 
Direct assimilation of coarse-scale observations is superior to a priori disaggregation. 
Independent assimilation of individual coarse-scale observations (method 3D-C1) can 
bring the overall mean analyzed field close to the truth, but does not necessarily improve 
estimates of the fine-scale structure. There is a clear benefit to simultaneously 
assimilating multiple coarse-scale observations (method 3D-Cm) even as the entire 
domain is observed, indicating that underlying spatial error correlations can be exploited 
to improve SWE estimates. Method 3D-Cm avoids artificial transitions at the coarse 
observation pixel boundaries and can reduce the RMSE by 60% when compared to the 
open loop in this study. 


Popular Summary: 

Satellites can indirectly measure the amount of snow that is present on the land surface. 
The measurement principle involves an analysis of the electromagnetic radiation in the 
microwave frequency range (centimeter wavelengths) that is naturally emitted by the land 
surface. Due to the characteristics of the satellite instruments and the measurement 
technique, the satellite observations can only provide information about the average snow 
amount within areas of about 25 km by 25 km. Such coarse-scale information is of 
limited value for water resources applications, particularly in mountainous areas. 
Complementary information on snow amounts can be obtained by using estimates of land 
surface characteristics (such as topography, vegetation, and soil information) along with 
precipitation (and other surface meteorological information) in a numerical model of land 
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surface processes. The model keeps track of the water and energy balance at the land 
surface and thereby also estimates snow amounts. Model estimates of snow can be 
obtained at relatively fine scales by setting up the numerical model on a 1 km by 1 km 
grid. 

Both the satellite observations and the model estimates are affected by necessary 
simplifications in the respective computational algorithms and by errors in the 
corresponding input data. An optimization technique known as “data assimilation” can 
be used to merge the information from the satellite observations and the land surface 
model. The resulting estimates are superior to the estimates from the satellite or the land 
surface model alone. By design, the data assimilation system distributes the coarse-scale 
satellite information onto the finer-scale model grid, thereby making the satellite 
observations more useful for water resources applications. 

In this paper, we investigate several options for implementing such a data assimilation 
system. One alternative is to disaggregate the coarse-scale satellite observations to the 
fine scale of the model prior to data assimilation or within the data assimilation system. 
Another alternative is for each model location to use only satellite observations that are 
local to that location or to also use observations from neighboring locations. We test 
these options for a 75 km by 100 km mountainous area in western Colorado, USA. We 
find that it is indeed possible to improve fine scale estimates of snow amounts through 
the assimilation of coarse-scale satellite observations. Our results indicate that it is best 
(i) to avoid disaggregating the satellite observations prior to data assimilation and (ii) to 
have satellite observations from the entire domain influence the assimilation estimate at a 
given location. The most sophisticated assimilation method avoids artificial boundaries 
around coarse-scale satellite observations and, for our study domain, reduces errors by 
60% when compared to the estimates provided by the land surface model alone. 
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